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An approach is proposed which allows to self-consistently calculate the structural and thermo- 
dynamic properties of highly charged aqueous colloidal suspensions. The method is based on the 
renormalized Jellium model with the background charge distribution related to the colloid-colloid 
correlation function. The theory is used to calculate the correlation functions and the effective col- 
loidal charges for suspension containing additional monovalent electrolyte. The predictions of the 
theory are in excellent agreement with the Monte Carlo simulations. 
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I. INTRODUCTION 

Over the span of the last century, colloidal suspensions 
have been the subject of intense theoretical and exper- 
imental study. The great effort is well justified by the 
importance that these systems play in industrial, biolog- 
ical, and medical applications. A practical problem that 
arises is how to stabilize suspensions against flocculation 
and precipitation, resulting from short range attractive 
van der Waals interaction. One approach is to synthesize 
colloidal particles with acidic groups on their surface. In 
aqueous environment these groups become ionized, re- 
sulting in repulsion between the macroions. 

Charge stabilized colloidal suspensions are an extreme 
example of a large asymmetry electrolyte. Both the 
charge and the size of the macroions are orders of mag- 
nitudes larger than those of other ionic species present 
inside the suspension. Typically a colloidal particle of 
radius 1000 A, will carry 10 3 - 10 4 ionizablc groups uni- 
formly distributed over its surface. The huge asymme- 
try between the macroions and the microions makes the 
theoretical investigation of colloidal suspensions a very 
difficult task [H, U H, 0| • The standard approach used 
to study these systems is based on the Primitive Model 
(PM), which treats solvent as a dielectric continuum of 
permittivity e. The interaction potential between the 
ionic species is taken to be composed of a long range 
Coulomb interaction and a short range hard-core repul- 
sion. Unfortunately due to the large charge and size 
asymmetry between the macroions and the microions, 
even for this simplified model the traditional methods 
of liquid state theory — such as the molecular dynam- 
ics simulations, the Monte Carlo (MC) simulations, and 
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the integral equations theories — prove to be only par- 
tially successful The huge number of counterions 
needed to ensure the bulk electroneutrality, allows the 
maximum charge asymmetry which can be studied using 
the present day computers to be around 100 : 1, while 
the important case of added electrolyte remains practi- 
cally unaccessible. Similarly, integral equation theories 
are plagued by convergence problems for strongly asym- 
metric electrolytes [5[- 

To obtain a more tractable description of these sys- 
tems it is, therefore, necessary to introduce further sim- 
plifications. This can be achieved by integrating out the 
microion degrees of freedom, leaving only a state depen- 
dent interaction potential between the colloidal particles. 
This defines the, so called, one component model (OCM). 
In spite of its apparent simplicity, the OCM requires 
knowledge of the effective macroion-macroion interac- 
tion, which implicitly depends on all the ionic species. 
Formally the potential can be obtained by explicitly trac- 
ing out the degrees of freedom of the microions of the 
PM @,@]. In practice, however, this coarse graining pro- 
cedure can only be accomplished by means of approxi- 
mate theories, such as the Poisson Boltzmann [y, 0) or 
the Ornstein-Zernike (OZ) equations, with appropriate 
closure relations @, SQ. Furthermore, to avoid compu- 
tational difficulties one usually assumes that the effective 
interaction potential is pairwise additive. This is quite 
reasonable at low macroion concentrations, however, care 
must be used when applying this assumption to more 
concentrated systems. As the concentration increases, 
the many-body correlations start to play an important 
role for both structural and thermodynamic properties. 
Assuming the OCM description with pairwise macroion 
interactions, there still remains a question of how to ob- 
tain the effective interaction potential. This has been the 
subject of many works [j], [T3, HH, EH- The difficulty in 
answering it is due to various factors, among which are 
strong correlations between the various particles and the 
huge asymmetry between the different ionic species - 
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forcing different approximations for different correlations. 
In the limit of large dilutions and small colloidal charge, 
a linearized Debye-Hiickel theory can be applied, and the 
pair potential takes a simple Yukawa-like form, known as 
the Derjaguin-Landau-Overbeek-Verwey (DLVO) poten- 
tial. For a system of colloidal particles of radius a, charge 
—Z D q, density p , and microions of valence Zi and bulk 
concentrations pi, i > 0, the DLVO potential is given by: 

/ z \ 2 P - K ( r ~ 2a ) 

P<r) =4 7—^ , (1) 

\{l + Ka)J r 

where lb = (3q 2 /e is the Bjerrum length, (3 — 1/fcsT, q 
is the elementary charge, e is the dielectric constant, and 

k = ^4:TT£ b J2i^ P' z i defines the inverse Debye screening 
length. From Eq. (1) one can see that, at this level of 
approximation, the role of the small ions is only to screen 
the electrostatic interaction between the macroions [l| . 

Even though the potential in Eq. (1) is restricted to 
low concentrations and small colloidal charges, the func- 
tional form of the DLVO potential can be extended to 
describe systems at moderate and high concentrations 
as well as large colloidal charge. To do this the struc- 
tural charge Z is replaced by an effective charge Z e ff, 
which accounts for the non-linear effects of counterion 
condensation [l], [H, [H, [lij]. In fact, it can be formally 
shown that non-linear short range correlations within 
the PM electrolytes can all be introduced into the DH 
theory by means of appropriate renormalization proce- 
dures 1 1 51 ] . The physical picture behind the charge renor- 
malization is that strong electrostatic attraction between 
the macroions and the counterions leads to their asso- 
ciation, so that from large distances (compared to the 
Debye length), a macroion can be viewed as carrying 
charge smaller than its structural bare charge. Both the 
macroion and its layer of condensed counterions can then 
be considered as forming a single entity of effective charge 
Z e ff. Once the non- linear correlations are taken into 
account through the charge renormalization, the DLVO 
pair potential, Eq. (1), can be used in the OCM de- 
scription to account for structural properties of colloidal 
suspensions. 

In this paper, we propose an ansatz which allows us to 
calculate both the thermodynamic and structural prop- 
erties of charge stabilized colloidal suspensions in a fully 
self-consistent way. This ansatz is based on a coupling of 
the renormalized Jellium model 0, [l7| with the OCM 
Ornstein-Zernike integral equations theory. From now 
on, we will only consider the case of aqueous monovalent 
electrolytes (zj = z± = ±1). 

II. THEORETICAL BACKGROUND 

Most of the theoretical work to obtain the effective 
charge of colloidal particles is based on the mean field 
Poisson-Boltzmann equation [13]. In many cases, the in- 
finite dilution limit is employed, and the problem reduces 



to that of a spherical macroion or an infinite planar wall 
immersed in electrolyte. For a more realistic situation 
of finite macroion concentration, the colloidal distribu- 
tion must be incorporated into the PB equation. To do 
this, one must solve the PB equation for a fixed macroion 
configuration from which the stress tensor and the force 
acting on each macroparticle can be calculated. Clearly, 
a numerical implementation of such procedure is very dif- 
ficult @, [H| • To have a more tractable approach, further 
simplifications are necessary. In this respect two approx- 
imations have proven to be particularly useful: the cell 
and the renormalized Jellium models. Before introduc- 
ing the new theory, we will make a brief review of the 
basic features of these approximations and discuss how 
the effective charges can be extracted from them. 



A. Renormalization models 

In the cell model, colloids are assumed to have a quasi 
solid-state like structure — macroions arranged in a form 
of a lattice. This allows us to consider one macroion in a 
corresponding Wigner-Seitz (WS) cell. A further approx- 
imation is to replace the polyhedral WS cell by a cell hav- 
ing the same symmetry as the macroion. The size of the 
cell is obtained from the overall macroion concentration. 
Because of the charge neutrality, the electric field must 
vanish on the surface of each cell, so that within this ap- 
proach there is no pair interaction between the colloidal 
particles. Nevertheless, the model is often used to calcu- 
late the effective macroion charges, which enters into the 
DLVO pair interaction potential. To obtain the effective 
charge, the non-linear PB equation is solved numerically 
inside the WS cell. The solution is then asymptotically 
matched to that of the linearized PB equation with an 
effective charge — the so called Alexander prescription 

The Jellium model captures the opposite limit in which 
the colloid-colloid correlation function is assumed to be 
completely disordered, g o{f) ~ 1 '9]. This approach is 
well suited for low density, weakly charged colloidal par- 
ticles. For strongly charged macroions, the Jellium ap- 
proximation fails to converge. Recently, Trizac and Levin 
have proposed a renormalization procedure designed to 
extend the validity of the Jellium approximation for 
strongly charged colloidal particles [l|, [4]. The renor- 
malized Jellium model relies on the concept of counte- 
rion condensation to determine the effective charge of the 
macroions. The method works as follows. One macroion 
with a charge Z Q is positioned at the origin of the coor- 
dinate system, the remaining macroions with their con- 
densed counterions are assumed to form a uniform neu- 
tralizing background in which the uncondensed counte- 
rions and coions move freely. Because it is not know 
how many counterions will condense onto the colloidal 
particles, the background charge density is not know a 
priori, but must be determined self-consistently. The dis- 
tribution of uncondensed counterions and coion around 
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the central macroion is assumed to be of the Boltzmann 
form, g j{r) = e^ ,3qz ^^ with Zj = ±1, where i/j(r) is the 
mean electrostatic potential around the central macroion. 
The electrostatic potential satisfies the modified Poisson- 
Boltzmann equation, 



VV(r) = | J2 PjZje-""^ ~ Z backPback (r) ) , 

'(2) 

where Z back and 

Pback 

(r) are the background charge 
and density, respectively. In the canonical ensemble — 
fixed number of all particles — pj 's are determined from 
the overall electroneutrality, while in the the semi-grand 
canonical ensemble (l7| . when the suspension is in con- 
tact with a salt reservoir at concentration c s , pj = c s . 
We note that Eq.Q would be exact if the electrostatic 
potential ip(r) on the right hand side of Eq.Q is replaced 
by the potential of mean force between the microion and 
colloid, w(r). In that case Zback would simply be the bare 
colloidal charge Z Q and Pback{r) = Po9oo{r), where p Q is 
the mean colloidal density. Unfortunately, there is no 
explicit way of calculating the potential of mean force. 
We are, thus, forced to identify w(r) « ip(r). This is 
permissible for monovalent ions in aqueous suspensions 
for which the electrostatic correlations between the mi- 
croions are small. The price for identifying w(r) m ip{r) 
is, however, a mandatory renormalization of the colloidal 
charge. Furthermore, one looses the direct identity be- 
tween the background density and the colloid-colloid cor- 
relation function. 

Within the renormalized Jellium approximation 
Pback(f) — p , and the bulk electroneutrality condition 
becomes 



■0qzjll) o 



ZbackPo — 0. 



(3) 



where V'oo is the Donnan potential which ensures the 
overall electroneutrality. In the canonical ensemble, we 
can take V'oo = 0. For a given set of parameters (includ- 
ing the background charge Zback), equation |(3J| can be 
solved numerically. Asymptotically, its solution has the 
form 



4>as(r) = Tpoo ~ 



Z efJ q e 



— K,(r — a) 



r(l 4- Ka) 



(4) 



In the semi-grand canonical ensemble [17| Pqipoo = 
— &rcsmh(Z e ffp /2c s ) and the inverse Debye length is 
k = y/%Tr£ b c s cosh^g^oo). In the canonical ensemble, 
ipoc = and k — \jAirtb(p- + p+). Eq. (0J allows us 
to calculate the effective charge Z e f / as a function of Z a 
and Zback- The self consistency condition is imposed by 
requiring that Z e ff = Zback, which determines the phys- 
ical value of the effective colloidal charge. It is important 
to note that unlike the cell model for which there is no 
pairwise interaction between the colloids, the macroion- 
macroion potential of the renormalized Jellium model is 
precisely of the DLVO form. 



To extend the renormalized Jellium model to larger 
concentrations, Castaheda-Priego et al. [2l[ have pro- 
posed to modify the uniform background density 
Pback(f) — Po, to account for the correlation hole around 
each macroions. These authors observed that for salt- 
free suspensions, simulations find that the colloid-colloid 
correlation function has the first maximum at r ss p~ x l 3 . 
They then suggested that this distance can be used to fix 
the size of the correlation hole between the macroions in 
salt- free suspensions [2l|. Castaheda-Priego et al. sug- 
gested that around each macroion there is an effective 
exclusion zone of radius r n = l/2p 1 ^ 3 , devoid of the back- 
ground charge. The factor of two is included in order to 
account for the fact that the exclusion zone is divided 
equally between the two macroions, see Fig. 1. The 
exclusion zone around each colloid is then taken into ac- 
count by replacing the usual uniform Jellium background 
density by a step function pback(r) — Po®(r — r h) in 
Eq. (2). Such procedure, however, still lacks the self- 
consistency, since the resulting effective charge can not 
be directly related with the correlation function, which 
is implicit in the form of Pback{r)- Furthermore, it is not 
clear how one can extend the above procedure to define 
the radius of the correlation hole for suspensions contain- 
ing additional 1:1 electrolyte. 



III. THE MODEL 

Although the renormalized Jellium model and its mod- 
ified versions allow us to calculate the effective charges, 
both theories lack in internal self consistency. In order to 
avoid this, it is necessary to find a way of calculate the 
effective charge and the correlation function g o(r) si- 
multaneously. To achieve this, we observe that the back- 
ground charge in Eq. ([2]) should be related in some way to 
the colloid-colloid correlation function. Unfortunately, as 
discussed above, within the renormalized Jellium model 
one can not identify the spatial variation of the back- 
ground charge directly with the g o{T). We note, how- 
ever, that g o{f) does carry the information about the 
size of the exclusion zone, which is approximately half 
the distance to the first peak of g 00 (r). In view of this 
observation, we will make the ansatz of identifying the 
background density variation with the rescaled colloid- 
colloid correlation function pback(r) = p goo(2r). This 
choice leads to a uniform background far from colloid 
pback{r) ~ poi while at the same time produces a corre- 
lation hole of appropriate size, Pback{r) ~ for r < rh. 

With this modification, the fully self-consistent Jellium 
equation becomes 



(f) = -47r4a 2 ( 



p+e 



-m - p _ e 4>(f)) +3r] Z back g 00 (2r), 

_(5) 

where we have defined the dimensionless quantities f = 
r/a, <(>(r) = /3eip(r), Z bac k = Z back £ b /a and rj = 
ina 3 p /3. As before, the effective charge is determined 
by the requirement that Z e ff = Zback- Eq. ([5]) is solved 
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FIG. 1: Two macroions of radius a separated by a distance 
r = 2x, at which the macroion-macroion correlation function 
has the first maximum. By symmetry, we can define around 
each macroion an effective exclusion zone of radius x. 

by an iterative procedure. We start with colloid-colloid 
pair correlation giV (r) = 1, and use this in Eq. ([5]) to ex- 
tract the corresponding effective and background charge 

■^e// = ^bacfe- The system is then considered in the OCM 
approach, with the interaction potential between the col- 
loids given by, 

1 / 7 \ 2 c-na(r — 2) 

I y (1 + KCL) I r 

where T = it/ a. In the semi-grand canonical ensemble, 
(no) 2 = (Kia) 2 (l + (Kres/ni) 4 ) 1 / 2 , («ia) 2 = 3rjZ eff and 
(K res a) 2 — 8ir£ c s a 2 . We then numerically solve the one 
component OZ integral equation with the Roger- Young 
(RY) closure to determine the new pair correlation func- 
tion g O0 •(f)- This function is then used as a new input in 
Eq. ([5|) to calculate the new effective charge Z\L. The 
procedure is iterated until the convergence is achieved, 

9oo{f) = goo (*")• In practice, only a few iterations are 
necessary to fulfill this condition. 

The Roger- Young closure is an interpolation between 
the hypernetted chain approximation (HNC) and the 
Percus-Yevick (PY) relation, with an adjustable param- 
eter a chosen so as to satisfy the thermodynamic self- 
consistency in the calculation of the isothermal compress- 
ibility 22]. In the salt-free case, the major contribution 
to the osmotic pressure comes from counterions, so that 
a is determined by imposing the requirement that 

where h OQ is the Fourier transform of the total correla- 
tion function. The first equality is the Kirkwood-Buff 
relation [23| . while the second one comes from approx- 
imating the microion pressure by the Jellium equation 
of state, PP = Z e ffp, and disregarding the weak depen- 
dence of the effective charge on the macroion density. 

For the case of large salt concentrations and moderate 
volume fractions — when the density dependence of the 
effective pair potential is weak — the pressure is given 
by that of the OCM [H,[25[, and the last equality in Eq. 



(J7J) is replaced by the OCM inverse compressibility, 

1 + PM0) = ^. (8) 

The OCM pressure Pocm can be calculated from the pair 
correlation function using the well known virial equation. 
It is important to note that in calculating the right hand 
side of Eq.®, the interaction potential must be kept 
constant [2J|. We use Eq.© to determine a in the RY 
closure when dealing with suspensions in contact with a 
salt reservoir at large concentration. 

In practice, due to finite discretization, the correlation 
function does not saturate at unity for large distances. 
Instead, it oscillates around 1 with a small amplitude. 
This creates difficulty for the solution of the PB equation. 
In order to ensure the correct long-distance behavior of 
ip(r), we set a cut-off distance f c , beyond which we force 
9oo{f) = 1. The value of f c is chosen such that \g o{f) — 
1| < 0.0025 for f> f c . 

IV. RESULTS 

In Fig. 2 we plot the macroion-macroion correlation 
functions calculated using the fully self-consistent Jel- 
lium (sc- Jellium) model developed above, and compare 
it with the results of the modified Jellium approximation 
of Castaneda-Priego et al. (m- Jellium) and with the MC 
simulations performed by Linse [26| for aqueous deion- 
ized suspension, K res = and colloidal volume 

fraction r; = 0.01. Three systems with coupling param- 
eters T = 0.3558 (a), 0.1779 (b), and 0.0445 (c), corre- 
sponding to particles of radius a w 160A, 40A, and 20A, 
respectively, were studied. All the calculations were per- 
formed in the saturation limit (very large bare charge). 
The corresponding effective charges are displayed in the 
graphs. As can be seen, both the sc- Jellium approach 
and m- Jellium show good agreement with the MC sim- 
ulations for T = 0.3558 and T = 0.1779, while for the 
lowest value T = 0.0445 the m- Jellium model seems to 
strongly overestimate the colloidal structure. 

In Fig. 3 we compare the correlation functions calcu- 
lated using the sc-Jellium with the MC simulations for 
various colloidal volume fractions at fixed coupling pa- 
rameter r = 0.3558 in the no-salt regime. Again, a good 
agreement with the MC simulations is found, for all the 
macroion concentrations. Surprisingly, this agreement 
seems to be better at higher volume fractions, dimin- 
ishing as the concentration becomes very low. Figure 
(4) shows the behavior of the effective charge as a func- 
tion of the colloidal volume fraction for the sc-Jellium 
(dashed curve), m- Jellium (dotted curve) and the origi- 
nal rcnormalized Jellium model of Trizac and Levin (solid 
curve). Although the qualitative behavior is the same for 
all three models, there is a significant quantitative vari- 
ation in the value of the effective charge. The effective 
charges predicted by the sc-Jellium lie between those of 
the m- Jellium and the renormalized Jellium models. 
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FIG. 2: Macroion-macroion correlation functions calculated 
using the sc-Jellium (dashed lines), m-Jellium (dotted lines) 
and MC simulations (solid lines) for a deionized colloidal sus- 
pensions at volume fraction rj = 0.01. The coupling parame- 
ters are (a) F = 0.3558 , (b) F = 0.1779 and (c) F = 0.0445. 



The real advantage of the sc-Jellium over the m-Jellium 
is that it allows us to accurately calculate the effec- 
tive charges and structures for suspensions containing 1:1 
electrolyte. At the moment this is the only theory which 
is capable of doing this for strongly charged colloidal par- 
ticles. The effects of non zero salt concentration on the 
macroion structure can be seen in Fig. (5), where the 
correlation functions for reservoir salt concentrations cor- 
responding to n res a = 1.0 and n res a =1.5 are displayed 
for various volume fractions. Unfortunately because of 
the difficulty of simulating these systems, no MC data is 



FIG. 3: Macroion-macroion correlation functions calculated 
using the sc-Jellium (dashed lines), and MC simulations (solid 
lines) for a deionized colloidal suspension with coupling pa- 
rameter F = 0.3558. From left to right, the volume fractions 
are given be r) = 0.08, r? = 0.04, rj = 0.02, r\ = 0.01, and 
77 = 0.005. 
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FIG. 4: Reduced effective charge as a function of volume frac- 
tion for a deionized colloidal suspension with F = 0.3558, as 
predicted by the renormalized Jellium model (solid curve) , m- 
Jellium (dotted curve) and the fully self-consistent approach 
developed in this paper (dashed curve). 



available. We see some very general trends for suspen- 
sions containing 1:1 electrolyte. As expected, increase 
of salt concentration leads to larger screening and loss 
of colloidal structure. In the salt dominated regime, the 
correlation functions become nearly independent of the 
macroion concentration. For these cases, both the effec- 
tive potential and the effective colloidal charge show very 
slow variation with the colloidal volume fraction — this 
also explains the weak variation of the correlation func- 
tions. Another remarkable feature is that at high salt 
concentrations, the colloidal structure is no longer im- 
portant for the computation of the effective charge, Fig. 
(6). The effective charge calculated using the the original 
renormalized Jellium model with a uniform background 
and the sc-Jellium are practically the same. 
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FIG. 6: Reduced effective charge as a function of volume 
fraction as predicted by: self-consistent approach with n r es = 
1.5 (dashed curve), and 1.0 (dotted curve) 
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FIG. 5: Macroion-macroion correlation functions calculated 
using the self-consistent approach, for suspension with T = 
0.3558 in contact with a salt reservoir at K res a = 1.0 in (a) 
and Kresd = 1.5 in (b). From left to right, colloidal volume 
fractions are 77 = 0.01, 77 = 0.02, 77 = 0.04, 77 = 0.01, 77 = 
0.005, 77 = 0.001 and 77 = 0.0005. As the volume fraction 
decreases, the correlation functions saturate, and the curves 
start to overlap. 



V. SUMMARY AND CONCLUSIONS 



We have developed a theory which allows us to 
sclf-consistently calculate both the thermodynamic and 
structural properties of aqueous colloidal suspensions 
containing 1:1 electrolyte. The theory is based on cou- 
pling the OZ equation with RY closure to the PB equa- 
tion with the renormalized Jellium approximation. For 
salt free suspensions, the predictions of the theory were 
compared to the MC simulations and were found to be 
in excellent agreement. The theory was then applied to 
study colloidal structure in suspensions containing 1:1 
electrolyte. Unfortunately, no simulational data is avail- 
able for these systems. Nevertheless, we expect that due 
to its internal self-consistency, the approach developed in 
this paper will remain very accurate for these systems as 
well. 
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